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FOREWORD 


The Industrial Mathematics Society is a professional organiza- 
tion whose objective it is to extend the understanding and application 
of mathematics in industry. Founded early in 1949, it marked the 
first organized attempt of any group to foster a closer relationship 
between mathematics and industry. 


It promotes the cause of mathematics in industry through 
monthly meetings at which technical papers are presented, through 
lectures by nationally prominent scientists and engineers, through 
panel groups which concentrate in specialized fields of knowledge, 
and through publication of worthwhile papers in its annual volume. 


Membership is open to engineers, scientists, and mathematicians 
who are concerned with the effective use of mathematics in industry. 
Annual dues are $5.00 per year including a subscription to the an- 
nual volume, Industrial Mathematics. Application forms and infor - 
mation may be obtained by writing to the Secretary, Industrial 
Mathematics Society, 100 Farnsworth, Detroit 2, Michigan. 


“Due to an editorial omission, there was no date printed on the 
outside cover of Volume 10, Part 1. In order to bring time and 
volume numbers back into proper phase, Volume 10, Part 2 has been 
labeled 1959-1960. This will permit the publication of Volume 11 
for 1961.” 


EXECUTIVE OFFICERS 
INDUSTRIAL MATHEMATICS SOCIETY 


President, Robert Davies, General Motors Corp., Research Labora- 
tories 

Vice President, Paul Loeber, Crysler Corp., Engineering Staff 

Secretary, Walter Kotoucek, Ford Motor Co., Engineering Staff 

Treasurer, Joseph Rackles, Ford Motor Co., Engineering Staff 
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Algebraic Fractions Applied to Logarithmic Functions 


By Douglas G. Anderson 
Associate with 

E. Kawecki, Design Engineer 
Franklin, Michigan 


ABSTRACT 


In this paper a method of interpolation is presented which per- 
mits abbreviation of six-place tables to the length of four-place 
tables without loss of accuracy. Tables adapted to the use of the 
method are appended. 

Also, a very brief formula is developed for use by a binary 
computer in direct computation of natural logarithms. The formula 
is correct within one unit in the tenth decimal place. 


INTRODUCTION 


Numerical computation by means of logarithms has dwindled 
since calculating machines have become common. Such machines, 
however, offer slight advantage for trigonometric computation, which 
requires frequent reference to tables and extraction of roots. It ap- 
pears that substantial savings of labor can be realized by shortening 
the tables. 

Such condensation requires an improved method of interpolation, 
since the principal reason for the bulk of the usual tables is the 
need to avoid serious loss of accuracy in interpolation. 


PLAN FOR THE TABLES 


To attain the greatest convenience for trigonometric computa- 
tion, each of the essential tables will appear on one double-page. 
The table of log N will give mantissas of the three-digit integers. 
The tables of log sin A and log tan A will have a tabular interval of 
ten minutes. An interval of .1 degree has much merit, but the 
number of entries required would crowd a page the size of those in 
this journal. The interval chosen is most convenient for work in 
degrees and minutes. 

l 
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The table of log sines will extend to 90°, and thus will contain 
the necessary log cosines. The table of log tangents will extend 
only to 45°, since log tan A = - log cot A = - log tan(90° - A). 

It is anticipated that a ten-inch slide-rule will be used for inter- 
polation. The formulas, however, will be sufficiently accurate to 
justify working to six places even where this is beyond the ability of 
the slide-rule. 

Linear interpolation will be used for angles larger than 30°, be- 
cause it is slightly easier, saves space in the tables, and is suffi- 
ciently accurate. 

The familiar functions, S and T, are given for angles less than 
2°, with the angles expressed in minutes. 


INTERPOLATION FOR LOGARITHMS OF NUMBERS 


Let N be a three-digit integer, and let n be a number less than 
one in absolute magnitude. It will be assumed throughout this dis- 
cussion that the base of logarithms is 10. 


2). 


log (N + n) = log N + log (1 +h 


(1) 


With the base 10, any number with the same sequence of digits has 
the same mantissa. The characteristic may be found by inspection, 
the mantissa of N may be read from the table, and the problem re- 
duces to finding the value of 


n 
d,, = log(l + 5) , (2) 


the interpolation function for numbers. 
From two well-known‘ series, both convergent when x?< 1, we 
obtain 


log(1 + x) = (log e)(x -% +3 -7 # s+) (3) 


and 





Therefore 











not 


we 














ALGEBRAIC FRACTIONS 3 


where the first term on the right is an approximation to d, and the 
error E, is 





_ Se. ae. ee 
E = (log e)(75 x -g* + 99% = eee), (6) 
Setting x, = - in eq. 5, and substituting in eq. 2, 
n loge 
d., “s 6 in +E... 


Since x, < = .01, E, < tf .000 000 04, completely negligible for a 
six-place table. 


INTERPOLATION FOR LOGARITHMS 
OF TRIGONOMETRIC FUNCTIONS 


Let A be a tabulated angle not less than 1°, and let a be an angle 
not larger than t 5’. 


sin(A + a) = sin Acosa + cosAsina 
= sin Acosa(l + cot Atana). 


Therefore 
log sin(A + a) = log sin A + log cos a + log (1+cot Atana). (7) 
Setting x. = cot A tan a, in eq. 5, gives 


tan a loge 





log(1 + cot A tana) = ah + 3 aa” Es. (8) 
For small angles expressed in minutes, tana = ror . Therefore 
eq. 8 may be written 
a log e 
= 22 a 9 
ds = 3438 tanA+2a* Bs: (9) 
Also 
teatA + a) = tanA+tana _ ns A (i + cot A tana 





1 - tanAtana 1 - tanAtana 


(tan A + cot A)tan + 
1 - tanAtana 





= tan A [li + 


tana 
sin AcosA - sin*A tana 





= tan A (1 + 


} 
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Therefore 


log tan(A + a) = log tan A + log(l + xt), (10) 


tan a 
sin A cos A - sin’ A tana‘ 





where x; = By eq. 5, where x = x;, 


x; log e 


log(1 + x,) ‘% tx + E, 


tan a loge 


~ gin A cos A +¢ (1 - 2sin? A)tana Ey 





. tan a loge 
3 Sin 2A + = tana cos 2A 





+ Ey 


ae | Seinie’ alll 
~ 1719 sin 2A +2 a cos 2A t 


a log e 
1719 sin 2A + 3a 





+E. + Ey. (11) 


The last step above was taken to separate the variables in the 
denominator of the interpolation function. The term thus introduced 
is, since log e = .4343, 


— ___.4343 a* (1 - cos 2A) — 
“ ~ 2(1719 sin 2A + 4 a)(1719 sin 2A + ¢ a cos 2A)’ 








which is closely approximated by 


4343 a*(1 - cos 2A) .4343 a? 





en ~- = =m came 2 
2(1719)? sin 2A (33a)? °° A- (12) 
Then eq. 10 may be approximated by 
4343 a 
log tan (A + a) = logtanA + (13) 





1719 sin 2A +3 a’ 


where E, and E.. are omitted. 
Thus the interpolation function for either log sine or log tangent 
may be represented by the same form, 


4343 a 


"Tay +4a (14) 

















—~— 











ALGEBRAIC FRACTIONS 5 
which represents the interpolation function for the log of the number 
(A+a). 
ACCURACY OF INTERPOLATION FUNCTION 
Eq. 7 may be written 


.4343 a 
3438 tan A + 5 a’ 





log sin (A + a) = log sinA + (15) 


where log cos a and E.~ are omitted. This omission introduces 
error. The sum of log cos a and E, has its greatest absolute 
value at a= - 5’. /E.! is greatest when !A! is smallest. Thus, at 
A=1° 10’, and a= -5’, x. = - .0714, and 


log cos a + E, = —.000 000 46 — .000 014 8 
= —.000 0153. 


When A = 2 10’, a similar computation yields -.000 002 6. 

Thus, if A > 2°, the error in log sine as computed by eq. 15 is 
less than t .000 0026. For large angles, the error is less than 
.000 000 5. 

Omission of E; and E, introduces error to eq. 13. Their sum is 
greater when a = 5’. When a = 5’ and A = 1°, x, = .0833, and 


E, + E. = .000 018 4 + .000 000 92 
= .000 0193. 


When A = 2°, a similar computation yields .000 003 4. For large 
angles, E, is completely negligible, but E. increases with A. At 30°, 
sec’ A = + and E. = .000 001 2. Thus, when 2° < A < 30°, the 
error in log tangent as computed by eq. 13 is less than ¢ .000 003 4. 

The errors computed above can be further reduced. E; and E. 
become less than .000 000 1 when A = 6°. When A < 6°, x, = x,=4, 
very nearly. Also, between 2° and 6°, d= .4343 . , approximately, 
and the major portion of E. or E; is .0362 (ay, Assuming a small 


adjustment, A, in the denominator of d, to cancel this error, 





aloge aioge _ loge (3) 
A-A A A 
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—_ 
from which A= 12 A° (16) 





The best choice for a in eq. 16 is that which minimizes the total 
error over the whole range of a. A study of the errors at A = 2° 
leads to 





_ 4343 a 
d, - + va (17) 
and 
_ 4343 a 
d, = R+ fa’ (18) 
in which 
Z = 3438 tanA - 000 510 (19) 
tan A 
and 
Zi .001 21 
R = 1719 sin 2A - sin 2A (20) 


Then, the error in log sine does not exceed t .000 001 when 
A> 2°. Also, the error in log tangent does not exceed + .000 001 2, 
from 2° to 30°. 


SLIDE-RULE ERROR 


A study of slide-rule accuracy indicates a probable error of .07% 
in a three-factor product as obtained by careful use of a ten-inch 
slide-rule. Accepting .10% as a conservative estimate, the probable 
error in evaluation of d is 


This error is in general much larger than all other errors. It is 
equivalent to an angular probable error of t .0010 a, which does not 
exceed .005 minute, or three-tenths of a second. This, in turn, is 
less than the probable angular error due to rounding off to six deci- 
mal places in the table of log sines for angles greater than 52°. 





anc 


on 





~ 


ALGEBRAIC FRACTIONS 
LINEAR INTERPOLATION 


For angles over 30°, A and R are omitted from the tables, and 
tabular differences are given. This saves space. Also, use of the 
linear interpolation function 


where D is the tabular difference, avoids the mental addition of 
{(A)+ 4a. The error does not exceed 


log sin 30° 5’ - 2 (log sin 30° 10’ + log sin 30°) = .000 001 8. 


SMALL ANGLES 


As A decreases below 2°, the interpolative errors rapidly in- 
crease. But, for small angles, sin A and tan A are very nearly pro- 
sin A oa © tan A 
A A 
are slowly varying functions of A, for which linear interpolation is 
accurate and easy. 

These functions have been frequently employed’, in the equa- 
tions 





portional to A. Therefore the functions S = 


log sin (A + a) = log (A + a) + S(A + a) 
and 


log tan (A + a) = log (A + a) + T(A + a). 


It has been most common to express (A + a) in seconds, but where 
the tabular interval in the main tables is ten minutes, it is better to 
express (A + a) in minutes. The special table appended to this paper 
gives S and T to six places without interpolation. 


INVERSE INTERPOLATION 


Solving the equation 
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for n gives the rule for inverse interpolation 


2! Se T-. 21 
4343 - 34, (21) 
Similarly, 
Zz d. 
a * 4343 - 24, (22) 
and 
Rd, 
a= 4343 - Fay (23) 


These formulas are easily memorized. 

Inverse interpolation for small angles requires that an approxi- 
mate value for the angle be known, in order to read S or T from the 
table. This approximation may readily be found on the slide-rule. 
Then 


log(A + a) = log sin(A + a) - S, (24) 
and 


log(A + a) = log tan(A + a) - T. (25) 


USE OF THE TABLES 


It is assumed that the user is familiar with logarithms. These 
directions will deal only with interpolation. 


1. Given a number, to find its logarithm: 
Disregarding the actual decimal point, think of the given num- 
ber as a number between 100 and 1000 with the given sequence 
of digits. From this number, subtract N, the nearest three- 
digit integer, to obtain n, a decimal fraction between - 0.5 
and 0.5. Read log N from the table, and compute 


.4343 n 


log (N + n) = logN + Wr Ta" 


The given number has the same mantissa as N + n. 
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2. Given a logarithm, to find the corresponding number: 
Find the nearest tabulated mantissa, log N, and read N. Sub- 
tract log N from log (N + n), the given mantissa. (The result- 
ing number, d,,, may be either positive or negative.) Then 


Nd, 
Nes? S——— 
4343 - 7 d., 
3. Given an angle between 2° and 30°, to find its log sine: 
From the given angle subtract A, the nearest tabulated angle, 
to obtain a, an angle between - 5’ and 5’. From the table, read 
log sin A, and Z. Compute the required log sine, 


.4343 a 


1 . 
Z+ 5a 





log sin(A + a) = log sinA + 


4. Given a log sine, to find the corresponding angle, between 2° 
and 30°: 
Find the nearest tabulated value, log sin A, and read A. Sub- 
tract log sin A from the given value, log sin(A + a), to ob- 
tain d.. Compute the required angle, 


Zd. 
ae °°* aa 


where a is expressed in minutes. 


5. Given an angle between 2° and 30°, to find its log tangent: 
From the given angle subtract A, the nearest tabulated angle, 
to obtain a, an angle between - 5’ and 5’. From the table, read 
log tan A, and R. Compute the required log tangent, 


.4343 a 
log tan(A + a) = logtanA + R+ H a 
6. Given a log tangent, to find the corresponding angle, between 2° 

and 30°: 
Find the nearest tabulated value, log tan A, and read A. Sub- 
tract log tan A from the given value, log tan(A + a), to ob- 
tain d,. Compute the required angle, 


R dad 


oe ee ae ee 


where a is expressed in minutes. 
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7. Given a small angle expressed in minutes, to find its log sine: 
Read S from the table of S and T, for the given angle, A. Com- 
pute log A by rule 1. Then 


log sinA = log A+S. 


8. Given a small angle expressed in minutes, to find its log tan- 
gent: 


Read T from the table of S and T, for the given angle, A. 
Compute log A by rule 1. Then 


log tan A = log A+T. 
9. Given the log sine of a small angle, to find the angle: 
Find an approximate value for the angle, A, by use of the slide- 


rule. With this value, find S from the table of S and T. Com- 
pute 


log A = logsinA -S. 
Find A, expressed in minutes, by rule 2. 


10. Given the log tangent of a small angle, to find the angle: 
Find an approximate value for the angle, A, by use of the slide- 
rule. With this value, find T from the table of S and T. Com- 
pute 


log A = logtanA - T. 
Find A, expressed in minutes, by rule 2. 


AN ACCURATE APPROXIMATION TO LOG, X 


Several methods have been developed for approximating log. x 
on high-speed computers. The more recent ones are modifications 
of a method due to I. J. Cherry’, which makes use of the relation 


loge x = (I + log, f)log.2, 


where x = 2'f, and} <f<1. This may also be written 
log.x = I log, 2 + log. f. 
In binary computers using “floating point” arithmetic, a number 


is expressed in terms of I and f. Thus the problem is to express 
log ..f, briefly and accurately. The series * 


' (1 + z) 
Be (i — 2) 





= 2z + (2/3)z* + (2/5)z° + (2/7)z7 +--+ (26) 





| 
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is a good starting point. Setting y = ; ., as y varies through 


the range from 7V2 to ¥2, z varies through the range from 
(3-2 Y2)to+(3 -2 2), which is smaller. Setting y= V2 f, we 


1 
obtain z = ip 3 7s and 
f+>y2 
log. x = (I - Zlog, 2 + log.y. (27) 


The series for log. y can be closely approximated by 


L = ( : —z)z = (a + c)z + (ab + cd)z* 


c 
1 -bz* 1-dz 
+ (ab* + cd?)z* 


+ (ab® + cd*)z7+---, (28) 


Initially, the constants may be chosen to make the first four 
terms of eq. 28 agree with the corresponding terms of eq. 26. Then 
a=1+ (1/18) ¥30, c= 1- (1/18) ¥30, b= 5 - (2/35) 730, and d= 3 


3 
7 
+ (2/35) 30. 
Comparison of higher terms of the two series shows that the 
error is 


logy - L = .01161 2° + .02578 z" + .0382z**° +--+. (29) 


When z = 3 - 2 V2 = .17157, the error is 15.98 x 107°. 

This error can be greatly reduced by small adjustments in a, b, 
cand d which leave the first three terms of the series in eq. 28 un- 
changed. Thus 


(a + Aa) + (c + Ac) =~artec. 


Oa = - Ac. (30) 

Also 

(a - Ac)(b + Ab) + (c + Ac)(d + Ad) = ab + cd. 

+ 

ae°-G = . ace =e 
And 

(a - Ac)(b + Ab)? + (c + Ac)(d + Ad)® = ab’ + cd’. 

acs . ME *+ ah m* tad + a@) Ad (32) 


(d + Ad)” - (b + Ab)? 
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Equating these two expressions for Ac, we obtain 


c(d - b) Ad 
b= ———————__—_—_—————__. (33) 
a(d - b) + (c + a) Ad 
These equations give exact values for all the changes, as fractions 
of the change in d, such that no error will occur in the first three 
terms of L. 


To determine the change in d, we expand the higher terms of L, 
neglecting powers of Ad higher than the first. 


AL = [(a - Ac)(b + Ab)* + (c + Ac)(d + Ad)® - (ab* + cd*)]z"” 
+ [(a - Ac)(b + Ab)* + (c + Ac)(d + Ad)* - (ab* + cd*)]z° 
Heee, (34) 


But, to a first-degree approximation, Ab = = Ad and Ac = - ac. Ad. 


Therefore, to a first-degree approximation, 


AL 
AD = (3(4? + b*) - 





— (d° - b*)]c 2” 


+ (ad? + b*) - 745 (at - vYJe at + - 


.2726 z" + 4673 2° + 5542z2%+--- 


The error equation is now 
E = log.y - (L + AL). 


In expanding this equation it is convenient to replace z by (3 - 2 V2)r 
= .17157 r. Then 


E = 43.76 x 107° | -272.6 Ad r* + (.3418 - 13.77 Ad) r® 
+ (.02234 - .480 Ad) r’ + ---]. 


Evaluating this for two values of Ad, say .0012 and .0010, and for 
several values of r, including r = 1, the most favorable value for Ad 
may be found graphically, as shown in the figure. When Ad = .001185 
the maximum absolute error is 1.09 x 107°. 

All that remains is to calculate the final values of the constants, 


and to put the expression in the most convenient form. This is 


o |» 
a|o 





log. y = i 7 FF 
. ; 
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Figure | 


The final values are 


. = 11 .245 464 629 1 = 933 140 917 51 
= 8.604 566 434 81 = 1,346 364 803 34. 
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COMMON LOGARITHMS. 
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The Present Theory of Switching and 
Some of its Future Trends 


By A. A. Mullin 
University of Illinois 


Our lot should not be so much one of wanting to construct ma- 
chines to do social-humanistic endeavors, but more one of wanting 
to keep man from doing any classical machine-like routines. 


0. Introduction: Since early time it has been recognized that the 
transportation of a note by a messenger usually entailed a long du- 
ration of time before the recipient received the note. Such a situ- 
ation placed nations with effective communication and transportation 
channels in a more favored political position than nations without 
such effective communication and transportation channels. The 
Appian Way and its connecting roads are a testament to the impor- 
tance of communication and transportation in the maintenance for 
400 years of a social-political-economic system called The Roman 
Empire. However, during the 19th century it became clear that the 
relative control of electricity, with its high velocity of propagation, 
offered a means of transporting messages to any place on earth ina 
matter of a fraction of a second. Such an idea soon led to the in- 
vention and development of telegraph, telephone, radio and tele- 
vision. In fact, this continuous emphasis on communication has led 
Wiener [1,2,3] and Shannon [4] to put forth the axiological proposi- 
tion that the two most important functions of each human are (1) to 
attempt more adequately to control his environment and (2) to com- 
municate with his fellow man. 

The creation of the dial telephone system and its complex 
switching apparatus necessary to join parties together, coupled with 
the speed and reliability of calculations that could be offered by 
computing machinery led to the development of a systematic pro- 
cedure for treating switching circuits. A switching circuit is an 
electric circuit composed only of componencs each having its elec- 
trical property described by being in either one or the other of two 
different states. This paper is an expository treatment of those 
mathematical aspects of switching theory that the author, at present, 
believes to be important. 
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1. History: The basic principles upon which the theory of switching 
circuits has been developed are those that treat the classical Aris- 
totelian syllogism [5]. During the 18th century the Swiss mathema- 
tician Leonard Euler [6] conceived the idea of treating syllogisms 
by means of logical diagrams and a few simple rules of manipu- 
lation. By a “logical diagram” is meant a diagram composed of dots 
and curves in which logical relations are denoted by such spatial 
relations that the necessary consequences of these logical relations 
are placed in clear evidence by the rules of manipulation. Near the 
middle of the 19th century the English mathematicians George Boole 
[7,8] and Augustus DeMorgan [9] conceived an algebra of conven- 
tional logic. This made syllogistic reasoning amendable to the more 
conventional algebraic manipulations. Toward the end of the 19th 
century the American mathematician and philosopher C. S. Peirce 
[10,11] considerably extended and developed Boole’s and DeMorgan’s 
algebra of logic. In addition, it is quite fair to consider Peirce as 
the inventor of the now familiar truth-table. In the last decade of 
the 19th century there appeared a book by the English logician John 
Venn [12] in which he considerably refined Euler’s idea of the use 
of logic diagrams. In fact, it was Venn who coined the term “Sym- 
bolic Logic.” In 1910, the Russian mathematician P. Erenfést [13], 
in his review of the translated copy of a book by the French mathe- 
matician L. Couturant [14] suggested the relation between two- 
valued logic and relay switching circuits. Supposedly, in 1935, the 
Russian mathematician V. I. Sestakov developed (but published in 
1941) Erénfést’s proposal. Further amplification of his results ap- 
pears in [15]. Between 1936 and 1938 there were numerous funda- 
mental papers published which developed the pure and applied as- 
pects of the theory of Boolean Algebra [16,17,18,19]. The Japanese 
engineers Nakasima and Hanzamus gave an application of Boolean 
algebra to switching circuits that preceeded the more thoroughly de- 
veloped works on the same subject by C. E. Shannon [20,21]. About 
the middle of World War 2, the classical paper by W. S. McCulloch 
and W. Pitts [22] was written which described an application of logic 
to a better understanding of the switching properties associated with 
the functioning of the human brain. In 1949, Shannon [23] published 
many new results concerning spring load distributions, contact load 
distributions and partially symmetric functions none of which are 
considered in [20,21]. In 1955, F. E. Hohn ana L. R. Schissler [24] 
published a paper on the use of matrices in the study of switching 
circuits. However, some of their results were already anticipated 
in [20]. In 1956, John von Neumann [25] published the results of his 
1952 lectures on Automata Theory that he gave at the California In- 
stitute of Technology. Certain parts of his paper are an extension 
and ramification of the results of [22]. In addition, he showed that 
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under certain conditions unreliable components could be coupled 
together so as to give arbitrarily reliable systems. At about this 
time D. A. Huffman [26], G. H. Mealy [27] and E. F. Moore [28] 
published their models for introducing the parameter of time into 
switching circuits. Some papers specifically directed at minimizing 
the number of components needed in the realization of Boolean func- 
tions are those of W. v.O. Quine [29,30], E. J. McCluskey, Jr. [31] 
and S. R. Petrick [32]. Moore and Shannon [33] consider the intui- 
tive ideas presented in [25] and apply them to relay contact net- 
works. In addition, they develop an interesting class of error- 
reducing network configurations called ‘hammock networks.’ D. E. 
Muller [34] applies the lattice theory of Garrett Birkhoff [35] to the 
analysis and synthesis of asynchronous switching circuits. A. A. 
Mullin extends the results of [33] to considerations with a parameter 
of time [36] and space [37]. Then he introduces the concept of a 
switching field [38] as a generalization of a switching circuit and 
metric reliability [39] as a unifying class of reliabilities. W. G. 
Kellner and he then treat a necessary condition for a minimum 
component of Boolean functions [40]. C. Y. Lee [41] considers a 
technique for representing switching circuits that is more in the 
form of a computer program. This procedure has the obvious side- 
advantage of giving a convenient means of specifying a measure of 
complexity for a Boolean function. At present, there are five mathe- 
matics or engineering textbooks of a general nature on the subject 
of logical design [42,43,44,45,59]. 


2. Classical Theory: 

2.1. Preliminaries: Section 2 of the paper treats the classical 
aspects of switching circuits. The word “classical” is used in the 
sense of involving the deterministic procedures of Boolean algebra 
for the study of switching circuits as proposed by Shannon [20,21,23]. 
In the paper we adopt the transmission function conventions [42]. 
Thus, ‘0’ denotes an open circuit, ‘1’ denotes a closed circuit, ‘+’ 
denotes parallel connections of contacts and ‘.’ denotes series con- 
nection of contacts. The classical theory of switching circuits is 
treated in two parts: 2.2. Classical combinational switching cir- 
cuits and 2.3. classical sequential switching circuits. For brevity, 
when in text, we drop the qualifying adjective ‘classical’ in referring 
to ‘classical combinational’ and ‘classical sequential’ switching 
circuits. ‘Classical’ is to be assumed present in remarks and only 
remarks concerning the ‘classical theory.’ By a combinational 
Switching circuit is meant a switching circuit whose electrical 
property is described, by means of Boolean algebra, only in terms 
of a sequence consisting of the states of each of the input variables. 
For example, a network composed exclusively of relay contacts is a 
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combinational circuit. By a sequential circuit is meant a switching 
circuit whose electrical property is described by means of lattice 
theory only in terms of the sequence consisting of the states of each 
of the input variables including the order in which each member of 
the set of input variables arrived at its state; for example, a net- 
work in which a contact for a relay also controls the excitation to 
the coil of the relay. 
2.2. Classical Combinational Switching Circuits: 

2.2.1. Preliminaries: The theory of section 2.2 is treated 
in two parts: 2.2.2 minimization theory and 2.2.3 unilateral iterative 
arrays. For practical reasons the treatment of combinational cir- 
cuits will not include contact-load and spring-load distribution 
problems [23], partially symmetric Boolean functions [31], bilateral 
iterative arrays, general multi-terminal contact networks [42], 
nor the general synthesis of non-series-parallel contact networks 
[24,31]. 

2.2.2. Minimization Theory: 

2.2.2.1. Preliminaries: By minimization theory we mean 
a procedure aimed at reducing the number of entities in a circuit 
that have a certain characterization. We consider only the problem 
of reducing the number of relay contacts or the number of electron 
tubes needed to realize a given Boolean Function. 

The theory of section 2.2.2 is treated in four parts: 
2.2.2.2. algebraic simplification, 2.2.2.3. symmetric functions, 
2.2.2.4. a residue test, and 2.2.2.5. the Quine-McCluskey reduction. 

2.2.2.2. Algebraic Simplification: Since each literal [29] 
in a Boolean function is to be realized by a component, one is moti- 
vated to apply identities for reducing the number of literals in a The 
representation for a given Boolean function. Eight such identities tact 
are the following: 
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Consider the circuit shown in Figure 1. Let T, denote the trans- 
mission function describing the electrical properties of the circuit 
shown. Then 


T, = A(C + BD) + (BC + B’)(A’ + D), or by (4) 
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Figure 1. A Combinational Circuit (9 contacts) 


= A(C + BD) + (B’ + C)(A’ + D), or by (1),(4) and (7) 
=C +AD+ A’B’. 


The simplified circuit for T, is shown in Figure 2. Thus, four con- 
tacts have been saved by some simple algebraic manipulations. 

If the original two-terminal network is non-series-parallel, the 
transmission function is found by tracing all paths between the two 
terminals. Consider the circuit shown in Figure 3. Then 


T. = wzx + wzw’y + wy’x’y + zx’y + zx’w’x + zy’zx 


or T. = z(x’y + x(w + y’). 


The simplified circuit for T, is shown in Figure 4. Thus, two con- 
tacts are saved by the algebraic manipulation. 

2.2.2.3. Symmetric Functions: An economical class of 
switching circuits is the one that realizes the class of symmetric 
functions. The basic concept of a symmetric function S of n varia- 
bles is that T = 1 if and only if some specified number of the varia- 
bles are 1. For example, S,, (x,, X,, X,) =X,X,+X,X,+X,xX,is a 
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Figure 2, Simplified Circuit for T, (5 contacts) 

















Figure 3. A Non-Series-Parallel Circuit 


symmetric function that takes on the value 1 if and only if at least 
two variables (i.e. either 2 or 3 variables) take on the value 1. Fig- 
ure 5 shows the basic circuit for generating all of the symmetric 
functions on three variables. By an iterative process the structure 
is easily extended to consider the case of any given number of vari- 
ables. 

2.2.2.4. A Residue Test: A residue test [40] is a con- 
venient means for determining the essential literals (29] in the 
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Figure 5. Generating all symmetric functions on three variables 


transmission function describing the electrical properties of a two- 

terminal combinational switching circuit [23], and thus indicates 

which literals are absolutely necessary in a realization of the 

switching circuit. 

The residue test is based upon a pair of identities of the 
Boolean algebra related by duality. These identities are 1 and 1’ 
below, where T is a transmission function. 





| (1) T = xx Ty + Xk Tox 
(1) T = [xx + To) + [Xk + Tx], 
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Consider the following three cases, in which the inclusion 
relation is the relation of possessing the same terms of standard 
' sum (in its canonic expansion) [42]: 
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Figure 6. A minimum contact circuit 


Case 1: Ty, = T,- Then T = x, Ty, + X% Tox = Tox = T,,- Here 
neither x; nor x’, is necessary for a realization of T. 


Case 2: 

(a) Ty. D> To, but Ti, A Ta (ie. Tu # Tox, and if To, = 1 then 
Tik= 1). Now To, + Ti, = To and To. + Ti, = Tix. Then T = xx (Tox 
+ Ty) +X Tok = Xk Tok + Tok. Here xi is not necessary for a re- 
alization of T. 

(b) To. > Ty but Ty # Tox (ie. Ty #T,., and if T,, =1 then 
T,x = 1). By reasoning similar to case (2a), the literal x, is not 
necessary in a realization of T. 


Case 3: None of the above. Here no reduction in the number of 
literals is possible. 
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If a circuit is found that uses each necessary literal 


(found according to the residue test) only once then the circuit is 
necessarily a minimum component circuit. Consider the circuit 
shown in Figure 6. The residue test shows that the only necessary 
literals are a, a’, b, b’, c, c’, d’. Since each of these occurs at most 
once the circuit shown in Figure 6 is a minimum contact circuit. 


2.2.2.5. Quine-McCluskey Reduction: The Quine- | 
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McCluskey [29,30,31] reduction procedure is essentially a system- | 


atic procedure for effecting all of the reduction identities of section 
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Figure 7. An Iterative Array 


2.2.2.2. Such a scheme is most useful to the circuit designer. How- 
ever, since no new mathematical ideas are concerned we will not 
consider the details. Such a detailed treatment is carried out in 
Caldwell’s book [42]. 

2.2.3.1. Unilateral Iterative Arrays: An iterative array 
is a collection of identical 2m-terminal contact networks coupled 
together as shown in Figure 7. Such an array is called unilateral to 
the right provided the circuit structure of the identical cells is such 
that a ground placed at the left of the entire array can propagate 
from each cell only to cells on its right. A typical transmission 
function that is economically realized by an iterative array is the 
following: T = 1 if and only if there is a single consecutive group of 
three operated relays. A typical cell for such a problem is shown 
in Figure 8. 

2.3. Classical Sequential Switching Circuits: 

2.3.1. Preliminaries: The theory of section 2.3 is treated 
in one part: 2.3.2 minimization theory. For practical reasons the 
treatment of sequential circuits will not include the recently de- 
veloped work by E. F. Moore [28] on sequential machines, the work 
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Figure 8. A Typical Cell 
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Figure 9. A “Lock-up” Circuit 
by D. A. Huffman [47] on delay polynomials, the work of C. C. Elgot 
[48,49] on automata theory, the work of W. G. Kellner [50] on iter- 
ative sequential circuits, pulsed sequential circuits [42], the recent 
work of Seiiti Huzino [51,52,53,54] on sequential machines, nor the 
work of J. P. Roth [55] on the application of algebraic topology to 
the synthesis of switching circuits. 
2.3.2. Minimization Theory: 

2.3.2.1. Preliminaries: We use minimization theory in 
the same sense as in section 2.2.2. The theory of section 2.3.2. i 
treated in one section 2.3.2.2. the Huffman model. For a minimiza- 
tion theory for other models see Moore [28] and Mealy [27]. We 
use the convention that a Boolean variable associated with a relay- 
coil is to be an upper-case letter and the Boolean variable associ- 
ated with contacts on that relay is to be denoted by the same lower- 
case letter. 

2.3.2.2. The Huffman Model: The Huffman model for a 
sequential circuit [26] consists of a number of combinational cir- 
cuits some of which control the excitation to relay coils (or other 
memory devices). However, some of the combinational circuits re- 
ceive a kind of ‘feedback’ from these coils in the sense that the coils 
actuate contacts in the combinational circuits. A simple example of 
a sequential circuit is the ‘lock-up’ circuit shown in Figure 9. The 
relays under our direct control at all time (input) are called ‘pri- 
mary relays’ and are denoted by an upper-case ‘x’ with the appro- 
priate subscript. All other relays are called ‘secondary relays’ and 
are denoted by an upper-case ‘y’ with the appropriate subscript. 
The output is denoted by a lower-case ‘z’ with the appropriate sub- 
script. Consider the following problem: Synthesize a two-input and 
one output relay sequential circuit whose output transmission func- 
tion is in the one state (i.e. grounded) if and only if the ordered-pair 
of inputs x,, x, form the connected sequence 00, 10, 11 and whose 
output transmission function is zero (i.e. ungrounded) for all other 
input-pair combinations. Double changes of input state are not per- 
mitted (i.e. the input transition 00 —-11 and 01 —410 are forbidden) 
in order to prevent ‘race’ conditions. 
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Figure 10, A Primitive Flow Table 


The classical synthesis procedure consists of the following: 

(1) Construct a primitive flow table with one stable-state per 
row. We know that at least one secondary relay must be used since 
there are only four possible input states. We denote stable-states 
by circled entries. An uncircled state seeks its correspondingly 
numbered circled state. Forbidden changes from a stable-state are 
denoted by a dash. Each stable-state determines a unique output - 
state and each stable-state depends only upon the previous stable 
state and the present input state x,,x,. Figure 10 shows a primitive 
flow table for the problem. 

(2) Next check the primitive flow table for inaccessible or 
equivalent states. Equivalent states at most differ in their notation 
in the sense that they lead to the same states for all input transi- 
tions. No states of the primitive flow table shown in Figure 10 are 
inaccessible nor is any pair of states equivalent. 

(3) In an effort to reduce the number of rows in the primitive 
flow table and hence possibly the number of secondary relays, a 
merger diagram is constructed. Such a diagram places in evidence 
all possible row mergers of the primitive flow table. Two or more 
rows can be merged, if column-wise, they have the same numerics 
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Figure 11. A Merger Diagram 


(disregarding whether an entry is circled or not) representing 
states. A node of the merger diagram corresponds to a row of the 
primitive flow table. The numbering of the nodes of the diagram 
corresponds to the stable state of the row of the primitive flow 
table. Join a straight line between each pair of nodes of the merger 
diagram that corresponds to rows that can be merged. A merger 
diagram for the problem is shown in Figure 11. 

(4) The merged flow table formed by making only the 4,5,6 
merger is shown in Figure 12. The dashed entries of the primitive 
flow table are chosen as the dotted numbers in the merged flow 
table indicate. 

(5) Now make an assignment of secondary states to the rows of 
the merged flow table (labeled a,b,c,d as in Figure 12) by means of 
a transition diagram. A transition diagram shows all interrow 
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Figure 12. A Merged Flow Table 
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A Transition Diagram 


transitions of the merged flow table by means of straight lines join- 
ing the nodes (which are identified with rows of the flow table) of the 
transition diagram. Figure 13 shows a transition diagram for the 
problem. Note that there is no difficulty with interrow transitions 
for the two left-hand columns of the merged flow table, since for 
these columns, there is no chance of going to any state but the de- 


sired state. 


No 


y, Y, — = O! 


0 100) 
01 100 | 


}+——__ + 





» 100] 10 | 


+ 


10 00 | lO} lO} 





EXCITATION MATRIX 


Figure | 


An Excitation and Output 
I 


Se 


EEE 


|| OF 


010] 


OUTPU T MATRIX 





Matrix 





A. A. MULLIN 


INPUT 





x 
1 Xo Y, - al 


x, v Y, 
+} 





a 
X, YH 
y, XY 





I 
rs 
i! 
z 


Figure 15, A Circuit for the Problem 


(6) Next write a secondary coil-excitation matrix and an output 
matrix each corresponding to the choice of secondary assignment to 
the merged flow table. A transmission function of 1 is associated 
only with state 3 of the merged flow table (this is a result of having 
been built into the primitive flow table). An excitation matrix and 
an Output matrix are shown in Figure 14. 


Hence, 
Y, = X, + X,Y, 
Y, = X,(,¥.+ X29 + Yi Ve) 
z= (x,x,y,). 


The underlined term in Y, is to remove the effect of a static 
hazard [46]. Otherwise, if the circuit was in stable-state 2 and we 
changed the input to x,,x, = 1,1, then the excitation to coil Y, could 
possibly be turned off momentarily. In general, the rule [46] is to 
include any prime implicants [29], necessary to prevent either a 
1-0-1 or 0-1-0 sequence of a transmission function controlling 
secondary relay coils. 

(7) Acircuit for the problem is shown in Figure 15. Note that 
we can use the minimization techniques of combinational circuits to 
get minimal circults for Y,, Y, and z. 


3. Stochastic Theory: 
3.1. Preliminaries: Section 3 of the paper treats the statistical 
aspects of the theory of switching. The word ‘stochastic’ is used in 
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the sense of involving the statistical procedures concerning fallible 
switching elements as proposed by von Neumann [25], Moore and 
Shannon [33], Mullin [36,37] and more recently [56]. 

The stochastic theory of switching may be divided into two parts. 
The first part is a calculus of probabilities built upon the super- 
structure of the classical Boolean algebra according to the usual 
postulates for probability theory, put forth by A. N. Kolmogorov [57]. 
Thus, if all probabilities are either zero or one then we return to 
the classical theory of switching circuits. The second part of the 
stochastic theory is a measure theory in a space of Boolean func- 
tions. These two problems are inherently different in nature. The 
first is suited to errors resulting in the m-triplication of Boolean 
variables and the second part is suited to errors resulting in the 
composition of different Boolean variables. The first part is appli- 
cable to relay contact networks and the second part is applicable to 
potential-level devices like electron-tubes, transistors and cryo- 
trons. In this paper we will consider only the first part dealing with 
a calculus of probabilities and called the ‘relay case.’ We will not 
treat the second part, but refer the reader to the works of von 
Neumann [25] and McCulloch [58]. The ‘relay case’ is divided into 
two parts: 3.2. relay stochastic combinational switching and 3.3. 
relay stochastic sequential switching. For brevity, when in text, we 
drop the qualifying adjectives ‘relay stochastic’ in referring to 
‘relay stochastic combination’ and ‘relay stochastic sequential’ 
switching. ‘Relay stochastic’ is to be assumed present in remarks 
and only remarks concerning the ‘stochastic theory.’ 

3.2. Relay Stochastic Combinational Switching: 

3.2.1. Preliminaries: The theory of section 3.2. is treated 
in four parts: 3.2.2. the model and reliability, 3.2.3. error-reducing 
configurations, 3.2.4. iteration processes, and 3.2.5. switching fields. 
The adjective ‘combinational’ is used in the same sense as for the 
classical case as described in section 2.1. 

3.2.2. The Model and Reliability: In the model for sto- 
chastic circuits each relay contact is defined by some binary infor - 
mation channel as shown in Figure 16. Here ‘a’ is the probability 
that the relay contact is open if the relay coil is deenergized. How- 
ever, for simplicity, we assume that a+b=1, that each and every 
relay contact is statistically specified by the same defining channel 
and that contacts function and malfunction statistically independent 
of one another. To this extent the probability p. of a contact mal- 
function is p, = l-a. 

By definition the reliability p of a contact or a circuit of 
contacts is expressed in terms of its probability of malfunction pe 
bp =-inp,. 
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Figure 16. A Defining Channel for a Relay Contact 


3.2.3. Ervror-Reducing Configurations: Let a and b be the 
parameters of a defining channel for contacts. When a relay con- 
tack is replaced by a two-terminal network of contacts on the same 
relay then the network of contacts has an equivalent defining channel 
as shown in Figure 17, where a is some probability function of a 
and b and § is some probability function of a and b (according to 
the conventions of section 3.2.2.). If 8< b and a <a then the net- 
work is said to be an error-reducing configuration. Figure 18 
shows two elements from the class of all error-reducing configu- 
rations. 

3.2.4. Iteration Processes: Let T‘*) denote a_ specified 
combinational circuit topology. Let Z be the set of positive inte- 


gers. Let T\J*'), j € Z, be the topology of the circuit formed from 
the topology T‘)), j « Z, by replacing each of its contacts with a cir- 


cuit of topology T) . Under these conditions T), d « Z is said to 
be the A-iterate of T''). The process under consideration is said to 
be the topological iteration of T’’’. It can be shown [33] that every 
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Figure 18. Two Error-Reducing Configurations 








iterate of an element from the class of error-reducing configu- 
rations is also an element of the class of error-reducing configu- 
rations. 

3.2.5. Switching Fields: Consider an isovolumetric topo- 
logical iteration process [37]. By this is meant the following: (i) an 
iteration of network topology using some element of the (non- 
vacuous) class of error-reducing configurations, and (ii) each net- 
work component has a reliability function defined on some non- 
vacuous index set of real-valued parameters such that the product 
of the total number of components by each of these parameters is a 
constant. 

By considering the circuit-reliability function, which is a 
real-valued function of network topology and component reliability, 
in the limit (if the limit exists) as each parameter of the domain of 
the circuit-reliability function approaches zero, one is led to an 
entity that may be called a homogeneous, isotropic combinational 
switching field. 

3.3. Relay Stochastic Sequential Switching: 

3.3.1. Preliminaries: The theory of section 3.3 is treated 
in two parts: 3.3.2. the model, and 3.3.3. stochastic matrices, error 
matrices and reliability matrices. 

3.3.2. The Model: By analogy to the case of classical se- 
quential circuits, the model for sequential circuits is to be a number 
of combinational circuits some of which control relay coils which 
have contacts in the combinational circuits. Then, by the techniques 
of section 3.2. we proceed to make each combinational circuit so 
reliable that the over-all sequential circuit has the desired relia- 
bility properties. We will not consider the case for which hazards 
are present. Such problems lead into the concept of a probabilistic 
prime implicant [36]. 

3.3.3. Stochastic matrices, error matrices and reliability 
matrices: Astochastic matrix is an exhaustive tabulation (analogous 
to the primitive flow table of classical theory) of each and every 
manner in which all the components function or malfunction, together 
with the transition probabilities between these states as determined 
by the defining channels for the combinational circuits (the defining 
channels for the components is assumed to be given). Let Ki; be 
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an entry of a stochastic matrix. Then ij = 1- H;; is the corre- 
sponding entry of the error matrix and Pj; = -1n (1-)H;; ) is the 
corresponding entry of the reliability matrix. If the entries of the 
reliability matrix are not too dispersed then P, = min Pj; is a use- 
ful measure of reliability for the sequential circuit. 


4. Future Trends: 
4.1. Preliminaries: Section 4 of the paper treats what the au- 
thor considers as future trends associated with the theory of switch- 
ing. The future trends in the theory of switching are treated in two 
parts: 4.2. trends in theory and 4.3. trends in applications. 
4.2. Trends in Theory: 

4.2.1. Preliminaries: The theory of section 4.2 is treated 
in four parts: 4.2.2. algebra, 4.2.3. mathematical logic, and 4.2.4. 
mathematical analysis and 4.2.5. topology. 

4.2.2. Algebra: Further work will continue on the applica- 
tion of the more recent results of lattice theory to the analysis and 
synthesis of asynchronous switching circuits. In this respect, the 
works of Muller [34] will be a most useful foundation. 

The use of Galois theory to study the general properties of 
delay polynomials [47] should produce some interesting results. 

4.2.3. Mathematical logic: New schemes of characterizing 
switching circuits, such as Lee’s [41] binary decision programs, 
will be sought. Recursive function theory should be useful in this 
respect (and incidently making a measure of complexity for Boolean 
function meaningful). More work will be done with the decision 
problem associated with automata. 

4.2.4. Mathematical Analysis: More work will be done to 
determine the properties of the class of error-reducing configu- 
rations. Mathematical tools will be developed (or rediscovered) so 
as to handle the case of potential-level stochastic switching circuits. 

4.2.5. Topology: Further applications of algebraic topology 
to the synthesis of switching circuits will be made. In addition, 
more use will be made of combinatorial topology in the synthesis of 
switching circuits. 

4.3. Trends in Applications: 

4.3.1. Preliminaries: The trends in applications are treated 
in three parts: 4.3.2. information transduction, 4.3.3. adaptive con- 
trol systems and 4.3.4. synthetic intelligence. 

4.3.2. Information Transduction: More use will be made of 
the basic properties of switching networks to provide error detec- 
tion and correction of transmitted messages. 

4.3.3. Adaptive control systems: Switching systems will 
play an important role in the analysis and synthesis of control sys- 
tems capable of adapting to themselves [60,61]. 
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4.3.4. Synthetic intelligence: Switching systems will play an 
important role when one is analyzing and synthesizing systems that 
show of the intuitive properties associated with intelligence [62,63]. 
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On Critical Shaft Speeds in Flexure of Spring 
Supported Automotive Drive Lines 


By H. K. Sachs 
Associate Professor 
Wayne State University 
College of Engineering 


The design of vehicular propeller shafts requires careful con- 
sideration of its dynamic behavior in flexure in addition to meeting 
manufacturer’s specifications. The need for a comprehensive study 
of critical shaft speeds led to a research contract awarded to Wayne 
State University’s College of Engineering by the Dana Corporation 
of Toledo, Ohio. The primary objective of this investigation is to 
determine the influence of elastic supports onthe shaft’s free vibra- 
tion frequency. 

The frequency study of rotating shafts is predicated on the fol- 
lowing assumptions: 


_ 


Shafts have a uniform and circular cross-section. 

2. The material is isotropic. 

3. Only flexural oscillations are considered, independent of 
torsional and longitudinal vibrations. 

4. The familiar relationship between shaft curvature and bend- 
ing moment (y” = - = is valid. 

5. The shaft is ideally balanced. 

6. Rotation of small shaft elements because of flexure and shear 
effects are neglected. 

7. The shaft is free of torque. 


The following notations are used throughout the article: 


< 
N 


Cartesian co-ordinates 

Circular frequency of shaft rotation 
Circular frequency of bounce motions 
Circular natural frequency of flexure 
Circular frequency of pitch motion 
Inertia moment of shaft cross-section 
Young’s Modulus 

Specific Weight 


were Qe 
i) 


~ 
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g Gravity constant 

U Shaft mass per unit length 

t Time 

L Shaft length 

a Argument of elastica per unit length 
(aL) = P Argument of elastica 

k,, k, Linear spring rate of supports 

K =k, +k, Resultant linear spring rate 

kb, , kb, Rate of flexural restraint at supports 


Kb = kb, + kb, Resultant rate of flexural restraint at supports 


= he Stiffness ratio of beam (Ed) and elastic sup- 
ports (k) 
Mx Bending moment at x 
Vx Shear force at x 
Qx Rate of change of shear force. (Also dynamic 
loading rate) 
Introduction 


The fundamental frequencies of vibrating beams having uniform 
and symmetrical cross-section may be obtained by equating the 
elastic restoring forces within the beam to the inertia force of mass 
at any arbitrary distance X from the support. The equations of 
motion of a shaft portion of unit length are: 


64 5? 
Lx Ee "2 a = 0 
5*z 6*z 
veh Get oe 


In Equations (1) the first terms denote the rate of change of the 
shear force Q, and the second term, the inertia load per unit length 
of shaft. A product solution of the form 


y = F, (x)o(t) z = F, (x)o(t) 
will satisfy the above partial differential equation. Since the shafts 
under consideration have a circular cross-section it will suffice to 


consider but one of Equations (1). 
A particular solution is known to be: 


n 
y= (C, sinh a,x + C,cosha,x + C, sina,x 
(2) 
+ C, cos a,,x)(sinw,t + cos w,t) 


(1) | 


_— 
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where C,, through C, are arbitrary constants to be determined from 
the conditions at the shaft supports. If we choose simple supports 
the following boundary conditions must be satisfied: 


2 2 
* +e -g SOY ‘ 
uw 9S vlan @ 8 re wae @ 0 bx? | x=1 = 0 (3) 


It can be shown that all arbitrary constants but C, must become 
zero in order to fulfill above requirements (Eq. 3) and hence the so- 
lution reduces to: 

* C, sin a, x (sin wt + cos w,t) (2a) 


and with a suitable choice of initial conditions this may be further 
reduced to: 


y, = C, sin a,x sin w,t (2b) 
for: 
Y,<-, = 0 aL =a, 2n---nga (4) 


where n is an integer denoting the n th noded mode of vibration. 
Substitution of Eq. (2b) in Eq. (1) yields the frequency equation: 


oe 
Wn  — 
and since 
a - aa 
L 
lt a IEnn 
n u L* 


for n= 1 we obtain the fundamental frequency: 


" = x ) IE (5) 

wn L? U j 
Expression (5) will assume the form of: 
/IE 

Wn = r. . where p77 (5a) 


for boundary conditions which are different from the assumed ones. 
Hence, the parameter p expressing the argument of the elastic line 
is dependent on the conditions of constraint at the shaft supports. 
The above expressions (5) are reasonably well known and yield 
nearly exact values of the fundamental shaft frequency (1,4). 
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Beams and shafts on elastic supports (2,3) 





In general, drive shafts rest on supports assumed to be rigid. 
If the shaft rests on very weak elastic supports, the natural fre- 
quency of shaft bounce motion can be determined from the relation: 


- /K 
@n -| m (6) 
whereby: m = uL = Shaft Mass 
For pitch motion we obtain: 
pa 
o, =V3V® (6a) 


Shaft assemblies which include universal or Hooke joints were found 
to flex between the joint yokes (the joint permits unrestrained bend- 
ing of shaft) if the shaft rigidity is considerably lower than that of 
the joints. However, vibrations of propeller shafts having a spring 
rate of approximately the same magnitude than that of the Hooke 
joints will cause flexural vibrations in the joints also. The exact 
determination of the effective shaft length “L” which exceeds that of 
the shaft proper is not always possible and will depend on the de- 
gree of rigidity of the supports. Hence, the elasticity of the support- 
ing structure must be expressed in the boundary equations. 
We shall consider two distinct boundary problems, namely: 


1. Shaft supported on springs of rate k, and k, (See Figure 1). 


2. Shaft systems according to (1) having also elastic con- 
straints of rate k,,, and k,,, at their ends (See Figure 2). 


( 
. 
q 


Figure 1, 
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It is clearly seen that the boundary conditions for Figure 1 re- 




















gid quire that: 
- wm. =o e2¥ = 0 
10n: x=0 “6x2 _ 
or 
6*y 
(6) M, 1 = 0 IE = ( 
Ox 2 
x=1 
also: (7) 
6° y : 
F eenit™ k, = -IE ra a = Vy-0 
5 
(6a) : - 6 7 
: y x=1 Xx k, IE 6x3 =e Vx 1 
und Upon substitution of the particular solution (Eq. 2) and its second 
- and third derivatives with respect to “X” in Eq. (7) we obtain four 
= homogeneous algebraic equations, in C,, C,, Cs, C4, viz.: 
‘ing 
oke | P se ¢&.+G . a) 
— * Rete a 
t of 
a a* (C, - C,) = 0 b) 
C, sinh aL + C, cos h aL - C, sin aL + C, cos aL = 0 c) (8) 
, C, (a8EI cos h aL - k, sin h aL) + C, (@ EI sin h aL - k, cos h aL) 
). 
on - -C, (a*EI cos aL + k, sin aL) + C, (a°EI sin aL-k, cos aL}=0 d) 


Cs a ae 








L 





Figure 2. 
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Since the characteristic values of the vibrating beam are inde. 
pendent of the initial conditions we may choose one arbitrary con- 
stant, say C, = 1. Writing the equations in dimensionless form, 
adopting the following notation: 


™ : EI 2 EI 
ie ee 


the remaining arbitrary constants C,, C;, C4 are found from Equa- 
tions (8b, 8c, 8d): 


C,=1 


pay [coth (coshp - cosp) - sinhp - sinp] + cosp 


he = 
*  p’p? [sinp cothp - cosp] - sinp 





_ p’p; [cothp(coshp - cosp) - sinhp - sinp] + cosp 
p’e2 [coshp - cotp sinhp] - sinhp 





, COSp - coshp 
sinhp 


Upon substitution of the expression of the arbitrary constants in 
Equation (8a) one obtains the frequency equation in terms of the 
parameter p (argument of the elastica) as a function of stiffness 
ratio p*, that is: 


pp, \cothp(coshp - cosp) - sinhp - sinp] + cosp 





p’p; [coshp - cotp sinhp] - sinhp 


pa, [cothp(coshp - cosp) - sinhp - sinp] + cosp (10) 
p*p? [sinp cothp - cosp] - sinp 





+ £O8P - coshp _ l 
sin hp pp” 


Next we shall consider the boundary conditions as shown in 
Figure 2 whereby, in addition to the support springs (k,, k,) two 
bending restraints (ky, > Ky) are attached at the shaft ends such as 


to prevent free angular displacement during flexure. The four 
boundary equations are: 





Afte: 


For 
with 


whe. 


in 
the 





CRITICAL SHAFT SPEEDS 5] 
_. &y| by" . oy 
(11) 
63 y 5?y by 
Xx nl - S «anise 
ae oe ox? | EE = ox] Kb: 








After differentiating the solution (Eq. 2) once, twice and three times 
and upon substitution in Eq. 11 we obtain again four homogeneous 
algebraic equations, viz.: 


k, C, +, 
EIC, -C;, 


a® + = 0 a) 


C, (a*EI cosh aL - k, sinh aL) + C, (a*°EI sinh aL - k, cosh aL) 
- C, (a® EI cos aL + k, sin aL) + C, (a°EI sin aL-k,cosaL)=0 b) 


(12) 
- odin & & Cc) 


C, (kp, cosh aL + Ela sinh aL) + C, (kp, sinh aL + ElacoshaL) 
+ Cs (kp,cos aL- ElasinaL) -C, (Kp, sin aL + ElacosaL)=0 d) 


Equations (12) reduce to a single equation in p, 


2, El 2_ El —_— El 
Ae S° kL Py EL Mb ELL 


For the purpose of brevity the final modal equation is presented 
without the intermediate steps. Then, 


2 
p' p? Ss DQ, ss GQ, a 2PPy Q, (13) 
1 "b, GQ, + DQ, 





where 


en p’ coshp - 1/p? sinhp _ 1/p§,, coshp + p sinhp 


1/p},, tanhp - p 
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B= p* cos p + 1/p? sinp 1/p},, COSP -psinp 














p*tanhp - 1/2 1/pi, tanhp + p 
p* sin p - (4) cnsp 1/p;,,, sinp + p cosp 
D = "Ta + = a 
p* tanhp - 1/p? 1/p;, tanhp + p 
p>? iP. ‘ 1 
E = sin h p (~x -—z) + coshp (p --773) 
Po, P, 4% 
p* p 1 
F = sinhp (—;—) -—sinp - p' cosp - =~ coshp 
P, p? Pr Pe 
4 2 1 ” 
p 
G = p*(sinhp + sinp) -p7 (coshp + cosp) 


As one can readily recognize Equation (13) is a transcendental alge- 
braic equation containing the parameter ‘“p” in an implicit form. 
“p” is found by an iterative method executed on a digital computer. 
Its value depends solely on the magnitude of the chosen stiffness 
ratios p* Af pi, p,, which, again, are functions of E, I, L and k,, 


K 2, k b, *b; 


DISCUSSION OF THE FREQUENCY EQUATIONS 





Both equation (10) and equation (13) yield infinitely many solu- 
tions for “p”. For a trial solution, the value of p = .1 is assumed. 
The following results are obtained for the first modes of a system 
as illustrated in Fig. 1. We obtain the following wave lengths: 


First mode p = 1.74 
Second mode p = 2.33 
Third mode p = 4,82 
Fourth mode p = 7.88 


It is not a coincidence that the difference between two successive 
values of p appear to approach 7 for modes above the third one. 
The hyperbolic terms “cos h” and “sin h” both tend towards infinity 
and “cot h” towards 1 for increasing values of P. Hence, their in- 
fluence on the frequency equation for the higher modes diminishes. 
The results of equations (10) and (13) are plotted in Fig. 3 and Fig. 4 
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Figure 3, 


respectively. Wave-lengths p corresponding to modes 1, 2, 3, and 4 
are plotted versus stiffness ratio p*. In Fig. 4 we show again wave- 
length p versus Pr, » the stiffness ratio with respect to the bending 
restraint, hence, we refer to Fig. 2. In order to plot the curves we 
assume 0° = .01 = constant. With this assumption the shear re- 
straint is given and only the bending restraint is variable. Of 
course once the wave-length p is determined we can substitute the 
value immediately into equation (5a) and obtain the circular fre- 
quency of flexure (See Fig. 9). Dividing the value by 27 we obtain 
the frequency in cycles per second. In the following we shall inter- 
pret the meaning of Fig. 3. 
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It is seen that for extremely small values of p, say less than .05 
the curve is rising steeply. Small values of the rigidity ratio imply 
very high stiffness values of the support springs. If we approach 
the limit which is zero we should obtain the frequencies which cor- 
respond to those of a simply supported beam. On the other hand if 
the rigidity ratio p* assumes an infinite value the condition implies 
supports of zero spring rates or, we may say, the supports are 
eliminated. It is seen that the first two modes expressed by wave- 
lengths p, and p, approach zero, that is they vanish. However, 
mode three and following approach a minimum value, e.g. mode 
three approximately 4.68 which corresponds to the wave-length of a 
“free-free” beam, or an unsupported beam. Let us consider the 
first two modes for weak elastic supports, that is, spring mounts 
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which result in a rigidity ratio of p* = .10. The first flexural mode 
degenerates into the state of bounce motion, the shaft being con- 
sidered as a rigid body. The second flexural mode degenerates into 
pitch motion. Yet, in both cases light flexure is thought to be super - 
imposed. over the motion pattern of the system as a rigid body. 
Hence, we recognize that the existence of external constraints in- 
creases the possible number of vibration modes by two, if we ex- 
clude motions in the direction of the length axis of the beam, and 
restrict ourselves to motion in a plane. 

It is advantageous to plot the deflection curves of the shafts for 
each particular mode. Thus, we obtain the configuration of the beam 
during vibration and we can immediately associate the mode to the 
number of node points obtained. The elastic lines for those modes 
which pertain to simple constraints (Fig.1) are illustrated in Fig. 5. 
Those pertaining to the combined shear and bending constraints are 
illustrated in Fig. 6. From these two illustrations we can immedi- 
ately conclude the effect of relaxation of constraints and observe 
that the two lowest modes of flexure degenerate, in effect to vibra- 
tion modes corresponding to bounce and pitch motion of the shaft as 
a rigid body. 

It remains to be seen what effect unsymmetrical constraints 
would have on the frequency of the system. Let us assume that the 
left hand constraints pb, and p,’, are unequal to the right hand con- 


straints Pre and p>. In general the frequencies observed tend more 


toward those values which correspond to the lower rigidity ratios 
regardless of the position of the weak elastic support with respect 
to the midpoint of the shaft. The elastic line of vibrating shafts 





















Figure 6, 


having unsymmetrical constraints is shown in Fig. 7, 7a, and Fig. 8. 
Figures 7 and 7a pertain to a shaft with unsymmetrical constraints, 
for example p; = .025, p} = .015 and p? = .015, p? = .025. Figure 8 








corresponds to unsymmetrical bending constraints, namely Pr, mus 
= 500, p. = 50 with constant shear constraints R* - g = 01. As 
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was expected the asymmetry of the support springs is reflected in 
the elastic line of the vibrating shafts. In all cases discussed, the 
condition of orthogonality of the natural modes, namely 
. 8. 1 
nts, J Yn (x) Vin (*) dx = 0 
re 8 


must hold. 
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SUMMARY OF RESULTS 























In conclusion of our discussion four important facts have been 
established: 


Re 


Elastic supports lower the natural frequency of beams sub- 
ject to flexural oscillation. 


The application of very weak support springs causes a de- 
generation of the two lowest flexure modes into vibration 
modes associated with bounce motion and pitch motion of a 
rigid mass. 


The effect of asymmetry of elastic supports is to reduce the 
natural frequency as to correspond rather with the weaker, 
than the more rigid support stiffness. Also, the deflection 
curve for any particular mode becomes asymmetrical. 

By varying the spring rate of right hand and left hand support 
between zero and infinity we can determine the natural fre- 
quencies of beams for any type of constraint which is con- 
templated including the following: 


a. Built-in at both ends. (p” = g* = pi = pi, = 0) 


b. Built-in at one end, simply supported at the other end 
%,, = p = pe = 0, Py, = 2) 
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c. Cantilever beam (%, =p' = Ph, = g? = 0) 


1 
d. Built-in, but movable at one end, free at the other 
2 a“ > im 2 = 2 _ 
= 9R =R "o,"7 1 
e. Built-in, but movable at one side, simply supported at the 


2 


other %,, = pe = 0, PL, = p; =o ) 
f. Built-in, but movable at both sides 


. 9a ae 
g- Built-in at one end, built-in and movable at other side 
,, = 9b, == 98 =) 
h. Simply supported at both ends (p{, = pi, =@ p? = = 0) 
i. Simply supported at one end, free at the other end 
= Po. =A =* » P, = 0) 


j. Free-free beam %,, = *, = p = p = ) 


Theoretical consideration and experimental evidence point to the 
desirability to suspend drive lines on elastic mountings. This will 
lower the critical vibration frequencies of flexure and will allow op- 
eration through these regions of instability in a virtually vibration- 
free manner. A forthcoming publication during 1961 of Soc. of Auto. 
Eng. will discuss the experimental aspects of this work in detail. 


ACKNOW LEDGMENTS 


The author is indebted to the Dana Corporation, Toledo, Ohio for 
their sponsorship of the drive line research. In particular he owes 
thanks to their Director of Research and Development, Mr. Philip 
Mazziotti, who made many helpful suggestions, and to Mr. Robert 
Warnke, also Dana Corporation, who undertook the programming 
work of the problem and execution of the same. Finally, thanks are 
expressed to all persons who wereconsulted during the research and 
contributed directly or indirectly in the compilation of the material. 


REFERENCES 


1, Den Hartog, Mechanical Vibrations, 4th Edition, McGraw-Hill, 
pp. 148, 432. 

2. Collatz, Eigenwertprobleme, Chelsea Publishing Company, pp. 
21, 42. 

3, Hohenemser - Prager, Dynamik der Stabwerke, 1933, 77. 

4. Timoshenko, S., Problems in Engineering, 3rd Edition, Van 
Nostrand, p. 324. 


























Applying a Matrix Partitioning Technique to 
Solving Steady State Vibration Problems * 


By James Duddles 
Pontiac Motor Division 


ABSTRACT 


This paper describes a matrix partitioning procedure which can 
be used for solving large systems of linear equations wnici have 2 
large quantity of zero terms. The use of a control matrix is intro- 
duced to save solution time for an example vibration problem. Fur- 
ther gains are experienced by rearrangement of the matrices as 
demonstrated on a sample problem. The techniques discussed have 
made possible solutions to problems which otherwise could not be 
handled on available computers. 


In vibration work, many systems can best be described by point 
masses connected together by weightless springs. A single mass in 
a system can move in six possible directions and therefore, re- 
quires six simultaneous equations to describe it. To simulate a 
complex system would require many mass points. For instance, a 
car body cannot be properly described with less than about 70 inter- 
connected point masses. Using the assumptions that the system is 
nearly linear and that steady state results are required, the car 
body would require 420 simultaneous linear equations to describe it. 

Digital computers lend themselves to solutions of linear equa- 
tions through the use of matrix operations. However, with limited 
storage capacity and the speed of operation of present computers, 
direct inversion of matrices of the order previously discussed are 
uneconomical. It is therefore necessary to use techniques which 
accomplish the solution through a method other than direct inver- 
sion. One solution which has been previously used is a partitioning 
technique’. Through the use of this method, advantage can be taken 
of the zero terms to cut down on the solution time. 

Partitioning is a process of dividing a large matrix into several 
smaller matrices. The process of solving a set of linear equations 


*I wish to acknowledge the assistance given by Mr, Fred Timpner 
of Pontiac Motor Division on both the original problem and on this 
paper, 
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Figure 1. An Example of Matrix Partitioning 


using the partitioning technique can be explained using the example 
equations shown in Figure 1. In the sample problem, the x’s repre- 
sent non-zero terms. The partitioning shown divided the large A 
matrix into square submatrices. However, this is not necessary 
for applying the partitioning technique. The equations in Figure 1 
can be represented by the three following equations where each term 
represents a submatrix. 


la.) A(l, 1) x ¥ (1) + A(1, 2) x ¥(2) + A(1, 3) x Y¥ (3) = B(1) 
1b.) A(2, 1) x ¥ (1) + A(2, 2) x ¥(2) + A(2, 3) x Y¥ (3) = B(2) 
lc.) A(3, 1) x ¥ (1) + A(3, 2) x ¥(2) + A(3, 3) x ¥ (3) = B(3) 


The solution for Y can be obtained in the following manner. 

The first restriction on the solution is that the matrix A(1, 1) 
must be of such a nature that its inverse can be found. On the as- 
sumption that an inverse, A(1, 1)™*, can be found for A(1, 1), the 
first step would be to multiply equation la.) by the inverse matrix 
A(1, 1)™* , to get equation 2a.). 


2a.) A(1, 1)™* x A(i, 1) x ¥(1) + A(1, 1)™* x A(1, 2) x ¥(2) 
+ A(1, 1)™* x A(1, 3) x ¥(3) = A(1, 1) x B(1) 
When equation 2a.) is multiplied by matrix A(2, 1), the result will be 
equation 2b.). Likewise, multiplying equations 2a.) by the matrix 
A(3, 1) will yield equation 2c.). 
2b.) A(2, 1) x ¥(1) + A(2, 1) x A(1, 1)™* x A(1, 2) x ¥(2) + A(2, 1) 
x A(1, 1)" x A(1, 3) x ¥(3) = A(2, 1) x A(1, 1)™ x B(1) 


2c.) A(3, 1) x ¥(1) + A(3, 1) x A(1, 1)™ x A(1, 2) x ¥(2) + A(3, 1) 
x A(1, 1)™ x A(1, 3) x ¥(3) = A(3, 1) x A(1, 1)” x B(1) 
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MATRIX PARTITIONING 


The three equations 2a.), 2b.) and 2c.) can be developed only if the 
matrices are of correct size for multiplication. This is the second 
restriction on the solution. Providing these restrictions have been 
met, equation 2b.) can be subtracted from equation 1b.) to give 
equation 3a.). In a similar manner, equation 2c.) can be subtracted 
from equation 1c.) to yield equation 3b.). 
3a.) [A(2, 2) - A(2, 1) x A(1, 1)™ x A(1, 2)] x ¥(2) 
+ [A(2, 3) - A(2, 1) x A(1, 1)™ x A(1, 3)} x ¥(3) 


= B(2) - A(2, 1) x A(1, 17° x B(1) 
3b.) [A(3, 2) - A(3, 1) x A(1, 1)™ x A(1, 2)] x ¥(2) 


+ [A(3, 3) - A(3, 1) x A(1, i)™ x A(1, 3)] x ¥(3) 
B(3) - A(3, 1) x A(1, 1)™ x B(1) 


The above equations can be simplified by the following substitutions: 


C(2, 2) = A(2, 2) - A(2, 1) x A(1, 1)™ x A(1, 2) 
C(2, 3) = A(2, 3) - A(2, 1) x A(1, 1)” x A(1, 3) 
C(3, 2) = A(3, 2) - A(3, 1) x A(1, 1)” x A(1, 2) 
C(3, 3) = A(3, 3) - A(3, 1) x A(1, 1) x A(1, 3) 
D(2) = B(2) - A(2, 1) x A(1, 1)™ x B(1) 
D(3) = B(3) - A(3, 1) x A(1, 1)” x B(1) 


" 


il} 


After making the above substitutions, the equations appear as fol- 
lows: 


3a.) C(2, 2) x ¥(2) + C(2, 3) x ¥(3) = D(2) 


3b.) C(3, 2) x ¥(2) + C(3, 3) x ¥(3) = D(3) 


Providing that the matrix C(2, 2) has an inverse, the equation 3a.) 
can be multiplied by the inverse matrix, C(2, 2) , to obtain equa- 
tion 4a.). 
4a.) C(2, 2)~* x C(2, 2) x ¥(2) + C(2, 2)™ x C(2, 3) x ¥(3) 
= C(2, 2)™ x D(2) 


When equation 4a.) is multiplied by the matrix C(3, 2), the result is 
equation 4b.). 
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4b.) C(3, 2) x ¥(2) + C(3, 2) x C(2, 2) x C(2, 3) x ¥(3) 
= C(3, 2) x C(2, 2) x D(2) 
Subtracting equation 4b.) from equation 3b.) will give the following 
equation 
5a.) [C(3, 3) - C(3, 2) x C(2, 2)™* x C(2, 3)] x ¥(3) 
= D(3) - C(3, 2) x C(2, 2)™ x D(2) 
When the following substitutions are made, the result will be equa- 
tion 5b.). 
E(3, 3) = C(3, 3) - C(3, 2) x C(2, 2)™ x C(2, 3) 
F(3) = D(3) - C(3, 2) x C(2, 2)~* x D(2) 


5b.) E(3, 3) x ¥(3) = F(3) 


From equation 5b.), the matrix Y(3) can be obtained. One method by 
which this can be accomplished is as follows. First, it is necessary 
to obtain the inverse of the matrix E(3, 3). Then, by multiplying the 
matrix F(3) by the inverse of matrix E(3, 3), the matrix Y(3) can be 
obtained. The matrix Y(2) can now be obtained from equation 6a.) 


by substituting in the matrix Y(3) and solving. Finally, the matrix 


Y(1) can be obtained by substituting matrices Y(2) and Y(3) in equa- 
tion la.). 

One important point needs emphasizing at this time. Notice that 
the C matrices can be placed over the A matrices without destroy- 
ing an A matrix that will be needed later on during this solution. 
Likewise, the D matrices can be placed over the B matrices, the E 
matrices over the C matrices, and the F matrices over the D 
matrices. This means that no additional space will be needed for 
the solution of the equations. 

There seems to be no advantage to this solution over an inver- 
sion solution until the submatrices containing all zero terms are 
considered. Although there are only 2 submatrices out of 9 which 
contain all zero terms, the number of matrix multiplications in the 
example problem drops from 17 to 11. However, it would be neces- 
sary to somehow know before the solution was started which sub- 
matrices were composed of all zero terms. One method of ac- 
complishing this would be to check all the submatrices for zero 
terms and set up a table. As an example, a table could be con- 
structed for the sample problem in Figure 1. The ones in the table 
will represent submatrices which do not have all zero terms and the 
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zeros will represent submatrices which have all zero terms. The 
table or control matrix for the A submatrices would appear as 
follows: 


























Figure 2. Control Matrix for A Submatrices 


By first looking in the control matrix, it is now possible to tell if a 
matrix multiplication is necessary. With the aid of the control ma- 
trix, it is now possible to cut down the time required for solution. 

In the application of partitioning to vibration problems, there 
can be considerable more savings. For instance, the control matrix 
previously described can perform an additional function. The zero 
terms in the control matrix indicate that no multiplications by the 
submatrix represented would be required. Therefore, there is no 
need to store the submatrix represented by the zero term. How- 
ever, the control matrix must now take on the function of locating 
the storage area of the non-zero submatrices. This can be accom- 
plished by storing the submatrices in a three dimensional array 
where the third subscript is used for location purposes. By storing 
the third subscript in the control matrix, the presence of a number 
greater than zero indicates that multiplication is necessary and tells 
where in storage to locate the submatrix. For example, the control 
matrix for the problem in Figure 1 might appear as follows: 


l 2 3 


Oc oe 0 
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Figure 3, Control Matrix Showing Storage Assignment 
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Using this arrangement, the storage required can be reduced for a 
problem, such as vibrations, containing some zero submatrices. 

At this point, a review of how the complete matrix is built up 
from the data for a vibration problem seems in order. The forces 
that a spring exerts cause two mass points to be coupled together 
and exert resistance to the movement of the mass points. If each 
mass is given the maximum of six degrees of freedom, the coupling 
terms will be contained in two six by six matrices which will be off 
the main diagonal. The resistance to mass movement terms will be 
contained in two six by six matrices, one for each mass point, which 
will lay along the diagonal of the matrix for the entire system. 
Since all the terms are contained in six by six matrices, the opti- 
mum size of submatrix for the system appears to be six by six. 
Using this size of submatrix allows storage of the fewest zero terms 
in the submatrices without getting the required control matrix too 
large and, also allows direct setup of the control matrix. There are 
two advantages of setting up the control matrix as the equations are 
developed. First, no searching is required later on to determine 
which submatrices are all zero and secondly, the complete matrix 
is never developed so the storage requirements for a given problem 
are lowered. 

When no storage is allowed for submatrices which contain all 
zero terms, there is a new problem introduced into the solution. 
For example, it is possible to end up with a control matrix similar 
to the one in Figure 4, where the numbers represent the non-zero 
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3) 3 x 4 x 12 
4 9 10 | x 
5 | 16 7 x 11 . x 
6 13 | x x 14 


























Figure 4, Sample Control Matrix 
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submatrices. During the process of solution, it will be found neces- 
sary to allocate storage for newly formed submatrices. The new 
submatrices required for the sample control matrix shown in Fig- 
ure 4 are indicated by the X’s. A general rule for determining in 
advance the additional storage required can be stated as follows. 
For non-zero submatrices above the diagonal, the submatrices in 
the column between the non-zero terms and the diagonal will be re- 
quired to complete the solution. For non-zero submatrices below 
the diagonal, the submatrices in the row between the non-zero terms 
and the diagonal terms will be required for the solution. Since these 
additional submatrices required for the solution cut down on the 
storage capacity and increase the solution time for a given problem, 
there is a need to keep them down to a minimum. 

One method of accomplishing a reduction in the number of ad- 
ditional submatrices required for the solution is through the use of 
row and column interchanges. Advantageous interchanges can be 
recognized easily when a matrix is being rearranged by hand. How- 
ever, letting a digital computer do the rearranging requires some 
organized method. In the following paragraphs, a description of the 
process used at Pontiac for selecting row and column interchanges 
will be described. 

The first step in arriving at the necessary interchanges is to 
make up a table showing the coupling terms connected to each point. 
Using Figure 4 as an example problem, the table would appear as 
shown in Figure 5. 


Point Coupling Terms 
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Figure 5. Coupling Table 


Using this table, a new order for the points can be developed. The 
new order is developed by taking the point number from the table 
and following the point number by its coupling terms. Before a 
number is added to the new table, a check is made to see if it al- 
ready appears in the table. If it does, the number is not added to 
the new table. The new table developed from Figure 5 would ap- 
pear as the first down column on left side of Figure 6. 
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First Second Third First Second Third Fourth 


Down Down Down Up Up Up Up 
1 1 1 1 3 6 6 
3 3 3 3 6 3 3 
5 5 5 6 1 1 1 
2 6 6 2 2 2 2 
6 2 2 5 5 5 +) 
4 4 4 4 4 4 4 

Figure 6. Rearranged Tables 


Starting with the first point, number 1, it must be placed in the new 
table. From the table in Figure 5, the coupling terms for point 1 
are found to be 3 and 5. Since neither of these terms are in the new 
table, they are added. The next point is 2 which does not appear in 
the new table and therefore, must be added. The coupling term for 
point 2 is 5 which already appears in the new table. The following 
point is 3 which already appears. Points 1 and 6 are the coupling 
terms for 3. Since point 1 already appears, the only addition to the 
table is point 6. The next point to be added to the table is point 4. 
No farther points will be added to the table since all the point num- 
bers already appear. The second down column on the left of Fig- 
ure 6 is developed using the same rules but instead of taking the 
points in numerical sequence, the points are taken in order from the 
first column in Figure 6. This process is repeated until the latest 
column matches the previous one. In the example problem, the third 
down column matches the second down column and this step of the 
process is terminated. 

The next step is essentially a continuation of the previous step. 
The major difference is that the new table is started at the bottom 
and proceeds to the top and the points are taken from the previous 
table starting at the bottom and working to the top. This process is 
continued until a new table matches a previous one. In the example, 
this procedure was stopped after the fourth up table as shown in 
Figure 6. 

The control matrix is then rearranged by column and row inter- 
changes to match the last table. The resulting control matrix would 
appear as shown in Figure 7. No additional submatrices are re- 
quired for the solution with this arrangement. However, before 
rearrangement, 8 additional submatrices were required for solution. 
Since the rearrangement was performed at the control matrix level, 
the actual quantity of numbers moved was small. In addition to re- 
arranging the control matrix, it is also necessary to rearrange the 
other matrices associated with the solution of the equations. For 
instance, the B matrix in the example in Figure 1 must also be 
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rearranged when the control matrix for the A matrix is rearranged. 
The resulting Y matrix will also be in a rearranged order and will 
require interchanges to put it in correct order. 

Through the use of partitioning with a control matrix and a re- 
arrangement technique, it is possible to accomplish the following 
improvements. First, the time required to solve a set of equations 
can be reduced by a large factor. Secondly, with a given storage on 
a computer, the number of simultaneous equations that can be solved 
will be greater. Finally, the combination of time savings and 
greater capacity makes possible the solution of problems which 
were previously impossible. 
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A Two Stage Inventory Model Applicable to a 
Provisioning Problem 


By Frank B. Quackenboss and Allen V. Butterworth 
Research Laboratories 

General Motors Corporation 

Warren, Michigan 


Arule is presented here under which two-stage inventory de- 
cisions can be made in the face of uncertainty. A common situation 
to which this rule applies is the one in which a buyer must deter- 
mine the correct amount to be purchased and stocked of some com- 
modity . 

The article to be stocked is assumed to be procurable in both 
fully -finished and partially-finished form. There are some obvious 
applications tothe stocking of spare service parts where the buyer’s 
decision is made at the end of a production run. 

In developing the rule, we assume that the decision maker has 
available to him 


1. A forecast. 

2. Reliable measures of the relevant costs of manufacturing and 
of obsolescence. 

3. Estimate the value of customer good will (returned or lost). 


The method of developing the rule to maximize expected profit 
follows the classical derivation of a single-stage inventory model 
with uncertainty. We extend the method, however, to two stages of 
inventory, i.e., partially finished and finished parts. 

Our forecast of demand must be explicitly in the form of the 
cumulative distribution function, F(n), which describes the proba- 
bility that the final outcome (sales) will be equal to or less than N. 
When the distribution function is absolutely continuous,’ the fre- 
quency function and the cumulative distribution function are re- 
lated by the expression 


. 
F(N) = J f(n)dn 


——— 


io ee , 
The assumption of an absolutely continuous distribution function 


is a good approximation for most economic situations. 
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Let: c = current cost of finished part per unit 
b = cost of partly finished part 
and 
a = cost to convert partially finished parts to finished pie 
parts at some future date (i.e., after the production 
run is terminated) 
R = selling price per unit 
M = good will of satisfied customer Fo: 
N, = optimum number of finished pieces stocked 
N,-N, = optimum number of partially finished pieces stocked. | 
Then N, and N, are given by j 
F(N,) = a+b-c 
. | for 
_(R+M) -a-b | sO 
F(N,) (R+M). a | sh 
| of 
provided the following inequalities are satisfied: th 
an 
R+M>a+b>c>b>Po | in 
and 


b 
(1 - ©) (M + R) 








hed 


ed. 
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The method of derivation is todetermine the change -in-expected- 
profit resulting from the purchase of an additional finished piece 
and to compare this change-in-expected profit with that resulting 
from the purchase of an additional partially finished piece. Each 
time we consider buying an additional piece, we select that action 
which gives the greater increase in expected profit. Thus additional 
finished pieces or partially finished pieces are added to stock until 
purchase of neither will increase profit. The change-in-expected- 
profit with the number of pieces purchases (both for partially - 
finished pieces and finished pieces) is linear in the cumulative dis- 
tribution of demand, F(N). 

The case of principal interest is shown in Figure 2 in which we 
plot the change-in-expected-profit vs. F(N) for both finished and 
partially -finished piece purchases. The change-in-expected-profit 
for added finished pieces is 


Pp 
oe = (R+M -c) - (R+ M) F(N,) (1) 
1 


and the change-in-expected-profit for added partially-finished 
pieces is 


aay = (R + M - a - b) - (R + M - a) F(N,) (2) 


For Figure (2), it has been assumed that 


R+M-.-croR?M-a-bd 
or 


a+boc (3) 


Inequality (3) means that the part which is bought now in finished 
form is cheaper than the partially finished part which is finished at 
some later date. If this inequality does not hold, all purchases 
should clearly be in partially finished pieces unless a = 0 (no cost 
of finishing) and b = c. In the latter case, there is no difference (for 
the purpose of this discussion) between purchase of finished parts 
and partial parts for later finishing. It has further been assumed 
in Figure 1 that 


(4) 
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Finished Pieces 
“~~ Slope - -(R + M) 


Partially Finished 
Pieces Slope = -(R + M- a) 








| F(N) 


F(N)) 





F(N.) 








Change of Expected Profit With Purchase of an Additional Piece < 
Lo ] 


Figure 2. 

Inequality (4) means that the total return from a sale (dollars 
plus good will) exceeds the cost of the item sold when partial parts 
are stocked for later finishing. If (4) did not hold, there would be no 
economic justification for stocking of partiaily finished parts. In- 
equality (5) means that there is some cost in converting partial 
parts to finished parts. 

Under the above assumptions, it can be seen from Figure 2 that 
purchase of finished pieces increases profit more than purchase of 
partial pieces until N, finished pieces are purchased. N, is de- 
termined from the equation 
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(R + M -c) - (R + M) F(N,) = (R+M-a -b) - (R+M - a) F(N,) 


or 


a+bD-c6 
a 








F(N,) = (6) 


| and the optimum quantity of finished pieces does not depend on 
either revenue from a sale or the value of customer good will. 
Quantities stocked in amount greater than N, should be partially 
finished pieces until additional stock no longer increases expected 
profit. As indicated in Figure 2, this occurs when the total quantity 
of pieces stocked (both as finished and partially finished pieces) 
| is N,, where N, is given by the equation 


(R + M - a - b) - (R + M - a) F(N,) = 0 





or 


R M-a-b 
3 Ie (7) 





F(N,) = 


Thus under the assumption resulting in Figure 1, N, finished pieces 
should be stocked. In addition, N, - N, partially finished pieces 
should be stocked. N, and N, are determined from equations (6) 
and (7). 

Even though the assumptions thus far expressed for Figure 2, 
namely 





R+Moart+bocr>bod 


are met, a case other than that shown in Figure 2 could arise. As- 
sume that “R + M” and “c” are the same as in Figure 2, but that “b” 
is increased and “a” is decreased so that “b” approaches ‘c” while 
“a+b” remains constant. Such a case is shown in Figure 3. As 
was found in Figure 2, Figure 3 shows that initially purchase of an 
ars additional finished piece increases expected profit more than pur- 
rts chase of an additional partial piece, but eventually a point is reached 
no | at which partially finished pieces are a “better buy” than finished 
In- pieces. 

‘ial | However, for the case shown in Figure 3, purchase of partially 
finished pieces is a “better buy” than purchase of finished pieces 
hat only in the sense that it causes profit to decrease less rapidly. At 
: of quantities for which the purchase of partially finished pieces is 
de - more desirable than purchase of finished pieces, stock of either ad- 
ditional finished pieces or additional partially finished pieces will 
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Partially Finished 
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Figure 3, 


decrease expected profit. Consequently, profit is maximized by 
stocking only finished pieces up to the quantity such that further 
purchases would decrease profit. Geometrically, this is the quantity 
corresponding to the value of the cumulative distribution, F(n), at 
which the ae vs F(n) line crosses the F(n) axis. Analytically, this 
is given by 


(R + M - cc) - (R + M) F(N,) = 0 
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If the inequalities 
R+Moa+t+bocobo 0 


hold, the cases shown in Figures 2 and 3 are the only ones possible. 
Case 1, the case of most practical interest, will arise if the 


r) 
finished part se vs F(n) line intersects the partially finished pieces 
a vs F(n) curve when = is positive. These two lines intersect at 


the point where 


a+b-c 
F(n) = ~ 
so that Case 1 will apply if 
(R+M-c)-(R+M) S&*P-<) 59 


or 


a(R + M -c) + (R+M)(c - a - b) = (R + M\(c - b) - ac > 0 


This reduces to 


(1 - —) (R + M) > a. 


om 


To complete the proof of the above statements, it is necessary 
only to show that the expressions for the change in profit with anad- 
ditional purchase be those which have been used above. These may 
be established in a straightforward fashion as follows: 


Let P(N,, n) = profit if N, finished pieces are purchased, given the 
demand is n. 
The expected profit, B(N,) is then given by 
- p00 
P(N,) = J) P(N, I, n) f(m)dn. 
0 


Now 


p 
Z 

} J 
" 


Rn - cN, 0O=n 


lA 
Z 


And 


P(N,! n) = (R - c)N - M(fn - N_) N<nas° 


1 
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Therefore 
ian Ni oo 
P(N, ) = f (Rn - cN,) f(n)dn+ { [(R - c)N, - M(n - N,)] f(n)dn 
N 
N, Lo) L<) 
= f Rnf(n)dn+ f RN, f(n)dn -cN,+M f (N, -n) f(n)dn 

‘ N; Ny 

dP(N,) 


_* RN, f(N,) - RN, f(N,) + SOR f(n)dn -c - M(N, -N,) f(N,) 
Ny 


1 





oo 
+M f{ f(n)dn 
N; 


= (R + M) J” t(n)an 
N, 


=(R+M -c) - (R + M) F(N,) 


Similarly, let 


} 
P(N,, N,!n) = profit if a total of N pieces are stocked, N, finished 


and (N,- N,) partially finished given the demand is n. 


Then 


P(N,,N,/n) = [Rn - cN, - b(N, - N,)]; O<n=N, 
= (Rn - cN, - b(N, - N,) - a(n -N,)]; N Sn=N, 
= (RN, - CN, - (a + b)(N, - N,) - M(n -N,)]; N,=n=* 


Jsing these three expressions for the expected profit resulting if N, 
finished pieces are purchased and partially finished pieces for 
(N, - N,) additional pieces are purchased is given by 


Ny 
B(N,, N,) = { [Rn - cN, - b(N, - N,)] f(n)dn 


Na 
+ { [Rn -cN, - b(N, -N,) - a(n -N,)] f(n)dn 
Ni 
«eo 
+ J [RN, -cN, - (a + b)(N, -N,) - M(m - N,)] f(n)dn 


Ng 


} 





A T' 





OSS 


f(N, 


} 
shed 
is n. 
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= f Rn f(n)dn + { RN, f(n)dn - cN, - b(N, - N,) 
N2 
Na 


- f a(n -N,) f(n)dn - f “a(N, - N,) f(n)dn 


N, N2 


-M f° (n-N,) £(n)dn 


Therefore, 


aT = RN, {(N,) - RN, {(N,) + f” R f(n)dn - b 
3 N 


2 


- a(N, - N,) £(N,) + a(N, - N,) £(N,) 


+ M(N, - N,) £(N,) + [~M f(n)dn 
Na 


=(R+M-a) f f(n)dn -b 
N, 
= (R+M -a) [1 - F(N,)] -b 


=(R+M-a-b) -(R+M -a) F(N,) 


aP 


Note that —— does not depend upon the value of N,, hence, the 


change in profit with additional purchase of partially finished pieces 
depend only upon the total of pieces purchased and not upon the way 
this total quantity is split between finished pieces and partially 
finished pieces. This fact made the graphical presentation given 


above possible. 


- {a f(n)dn 


Na 


er 


' 
: 
: 
f 














INDUSTRIAL MATHEMATICS is published semi-annually by the In- 
dustrial Mathematics Society. The advance subscription price is 
$4.00 per annual volume of two issues. Back issues (one volume 
per year) are available at $3.00 per copy. Orders accompanied by 
check or money order should be addressed to Industrial Mathe- 
matics Society, 100 Farnsworth Street, Detroit 2, Michigan. 


INDUSTRIAL MATHEMATICS contains originalarticles dealing with 
the application of mathematics in industry. Persons having such 
manuscripts for publication should submit them in duplicate to the 
Editor, INDUSTRIAL MATHEMATICS, 100 Farnsworth Street, 
Detroit 2, Michigan. Any questions on editorial matters should be 
similarly addressed. 


Manuscripts should be typed double-spaced on one side of 85 by 
11 inch paper in complete and finished form precisely as they are to 
appear, and sent — packed flat — to the Editor. Where possible, 
Figures should be numbered consecutively in arabic numerals. 
Tables should be numbered in arabic numerals, and both figures and 
tables should have captions. 


Authors will receive galley proofs for correction purposes. One- 
hundred reprints without covers will be furnished free. 


Author’s institutions or companies are asked to pay a publication 
charge of $7.50 per printed page for articles accepted for publi- 
cation, 











‘ 
; P 
: 
* 
$ ca 
; s 
. > A; 
. - 
; ° 
s 
. 9 
‘ 4 = ; 
¢ ‘ : : 
. . 
. % is ; c 
a e 
a 
: j 
E * 
{ 
. 
oe : 
Cee \ 
a 
. be 
| F 
> 
. 5 . 
’ 
* FS : ’ 
Z . 
. 
. 
; 7 
i . 
e 
. ’ 
_—, 
. 
. 
. 
. . 
' - 





